Self-consistent theory of molecular switching 
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We study the model of a molecular switch comprised of a molecule with a soft vibrational degree 
of freedom coupled to metallic leads. In the presence of strong electron-ion interaction, different 
charge states of the molecule correspond to substantially different ionic configurations, which can 
lead to very slow switching between energetically close configurations (Franck-Condon blockade). 
Application of transport voltage, however, can drive the molecule far out of thermal equilibrium 
and thus dramatically accelerate the switching. Fhe tunneling electrons play the role of a heat bath 
with an effective temperature dependent on the applied transport voltage. Including the transport- 
induced "heating" selfconsistently, we determine the stationary current-voltage characteristics of the 
device, and the switching dynamics for symmetric and asymmetric devices. We also study the effects 
of an extra dissipative environment and demonstrate that it can lead to enhanced non-linearities in 
the transport properties of the device and dramatically suppress the switching dynamics. 
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I. INTRODUCTION 

The apparent limitations of the silicon-based technol- 
ogy on the way to further acceleration and miniatur- 
ization have prompted active research into alternative 
electronic architectures^ In particular, molecular elec- 
tronics holds a lot of promise because each molecule, be- 
ing only about a nanometer in size can in principle per- 
form such non-trivial operations as information storage* 
or electrical current rectification^ Since molecules, even 
intricate ones, can be mass produced by means of well- 
controlled chemical synthesis, one expects them to be 
less susceptible to the issues of disorder that plague the 
silicon-based electronics below the 10 nm scale. The 
ultra-miniaturization that molecular electronics affords, 
however, also leads to the problem of connecting the 
molecular elements among each other, as well as of the 
necessary interfacing with large-scale conventional elec- 
tronics. Indeed, early on, this problem has caused many 
difficulties in reproducing results from one device to 
another. 5 However, recent advances in fabrication^ as 
well as better theoretical understanding of physics and 
chemistry at the point of contact demonstrate that this 
difficulty is not fundamental and promise to make reliable 
and reproducible molecular junctions a reality. 

There is, however, a fundamental difference that dis- 
tinguishes the molecular devices from the conventional 
semiconductor ones. For a molecule to perform its unique 
function, it has to be well isolated from most environmen- 
tal influences, except for the (metallic or semiconducting) 



contacts that are required to access it. Under standard 
operation of the device, the chemical potentials differ by 
the value of the applied transport voltage V mutiplicd 
by the electron charge e, and thus the environment that 
the molecule experiences can not be considered as equi- 
librium if the voltage is greater than the temperature, 
eV > ksT (fee being Boltzmann constant). Therefore, 
to determine the behavior of a molecular device under 
such conditions, one needs to determine selfconsistently 
the influence, e.g. of electrical current on the molecular 
dynamics, and vice versa, the influence of non-thermal 
vibrations or electronic excitations of the molecule on 
the current. This is very different from the conventional 
electronics where devices are rarely driven out of ther- 
mal equilibrium far enough to significantly affect the per- 
formance (exceptions are the non-linear devices, such as 
Gunn diod). 

One of the most promising and interesting molecular 
devices is a switch, which can be used for information 
storage. Switching has been observed experimentally in 
several molecular junctions. 8 Proposed theoretical expla- 
nations for switching range from (a) large and small- 
scale molecular conformational changes, (b) changes in 
the charge state of the molecule, or (c) combination of 
the two, or "polaronic" . The purely electronic switching 
mechanism (b), while possible, appears quite impracti- 
cal since it would require a separate contact in order 
to change the charge state of the part of the molecule 
that would play the role analogous to the floating gate 
in flash memory by electrostatically affecting the "chan- 
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nel" current. The switching mechanisms (a) and (c) 
upon closer inspection turn out to be fundamentally the 
same, since in order to be able to switch and read out 
the conformational state of the molecule electronically 
there necessarily has to be a coupling between the elec- 
tronic and ionic degrees of freedom. The dynamical sta- 
bility of the "on" and "off' states in these mechanisms 
is achieved due to the collective nature of the states, 
which now involve not only the electronic occupancy but 
also all the positions of the ions in the molecule. Thus 
the change of the charge state of the molecules is ac- 
companied by the ionic rearrangement, which for strong 
enough electron-ion coupling can dramatically slow down 
the charge state switching. This is the essence of the 
Franck-Condon "blockade." 9 ' 10 ! 11 ' 12 By chemically engi- 
neering molecules with strong electron-ion coupling and 
soft (low frequency) vibrational modes one can achieve 
arbitrarily slow equilibrium switching rates. 

In order for molecular memory element to be useful, 
it has to have a long retention time (slow switching rate 
in the absence of any drive) but fast write time, i.e. it 
should be possible to accelerate the switching rate by gate 
or transport voltages. It is easy to see that in molec- 
ular switches with polaronic mechanism the transport- 
driven switching acceleration occurs naturally. As soon 
as the transport voltage exceeds the vibrational energy 
quantum, eV > Hluq, additional transport, as well as 
switching, channels open, which correspond to electron 
tunneling on and off the molecule with simultaneous ex- 
citation of vibrational quanta i 10 ' 12 Moreover, enhanced 
charge fluctuations on the molecule effectively "heat up" 
the molecule, further increasing the current through the 
device. This leads to a positive feed-back loop which sat- 
urates when the energy transferred to the molecule from 
non-equilibrium tunneling electrons exactly balances the 
energy transferred back from the molecule to electrons. 
As a result, the stationary switching rate can vastly ex- 
ceed the equilibrium switching strongly suppressed due 
to the Franck-Condon physics. 

Most of the molecular devices studied experi- 
mentally so far have been weakly coupled to the 
leadsi 13 ' 14 ' 15 ' 16 ' 17 ' 18 ' 19 This corresponds to the bare tun- 
nel broadening TiT of molecular electronic levels smaller 
that the energy required to excite one oscillator quan- 
tum (phonon) tko . The single-electron effects play a 
crucial role in this case. They are well theoretically de- 
scribed by a model of a single-electron tunneling (SET) 
device coupled to a single-mode harmonic oscillator, de- 
veloped mostly in the context of nanoelectromechani- 
cal systems. In the strong-coupling regime, when the 
electron-ion interaction energy E p (defined below) ex- 
ceeds Tiu)q 1 the physics is governed by the Franck-Condon 
effect, i.e. when the tunneling of an electron onto the 
molecule with the simultaneous emission or absorption of 
several phonons is more probable than elastic tunneling. 
The current as the function of voltage exhibits steps sep- 
arated by foaJo/e ) 9 ' 20 ' 21 ' 22 and the non-equilibrium elec- 
tronic heating of the molecular vibrational mode leads 



to self-similar avalanche dynamics of current with the in- 
tervals of large current alternating with the periods of 
strongly suppressed current pi - 

In this paper, we study the case of "slow" phonons 
at strong coupling, T ^> lu for eV > ^qi 11 ' 12 ' 23 ' 24 ' 25 
The physical distinction between this case and the one of 
"fast" phonons, r <C ^o, can be understood in the follow- 
ing way. For fast phonons, every electron tunneling event 
occurs over many oscillator periods. Thus effectively elec- 
trons can only couple to (or "measure") the energy (i.e. 
occupation number) of the oscillator. 26 ' 27 In the opposite 
regime, T <C ojq, electron tunneling is fast, and thus elec- 
trons are sensitive to the position of the oscillator. There- 
fore, in the former result of electron tunneling, 
the oscillator density matrix becomes close to diagonal in 
occupation number basis (and thus non- classical), and 
in the latter case, it is nearly dia gon al in the position 
basis (and thus classical). In Ref. [l2| it has been rigor- 
ously demonstrated for arbitrary coupling that the con- 
dition for the onset of the classical (Langevin) dynamics 
is given by min(?ir, eV) 3> Tiluq. Even at weak coupling, 
E p < hujQ, if a high enough bias is applied between the 
leads, the oscillator dynamics becomes non-trivial, with 
the possibility of switching between stationary states of 
different amplitudes^ At strong couplings, E p > Tiwo, 
there is another kind of multistability that appears at rel- 
atively small voltages, eV < E p - the system can switch 
between the states corresponding to (approximately) 1 
and electrons on the molecule. This multistability and 
switching can be described within the generalization of 
the Born-Oppenhemier approach to open systems.— In 
the metallic case the appearance of the multistability and 
the current suppression as a function of the bias voltage 
is associated with a discontinuity of the current (when 
the cotunnelling is neglected)^ 

The slow (or "classical") phonon strong coupling case 
is attractive since besides switching between the differ- 
ent charge-ion states, it allows a read-out of the state 
by means of cotunneling transport through the molecule. 
In cotunneling, the charge state of the molecule changes 
only virtually for a period of time determined by the 
energy uncertainty principle. This time can be much 
shorter than the vibration period, and thus the ionic 
configuration and the average charge occupancy need 
not change. On the other hand, in sequential tunneling, 
the tunneling events between the leads and the molecule 
are energy-conserving, with the rates determined by the 
Fermi's Golden rule. Typically, cotunneling currents are 
much smaller than sequential ones since they are higher 
order in the tunneling matrix element. However, if the se- 
quential tunneling is strongly suppressed by the Franck- 
Condon physics, the cotunneling, which needs not be af- 
fected by it, may dominate. In the case T < loq and 
strong electron-ion coupling the role of cotunneling was 
recently studied in Ref. [29|, where it was found that while 
it does not destroy the Franck-Condon blockade, it can 
dramatically affect the low-voltage current and current 
noise, as well as the vibrational dynamics. 
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The purpose of this work is to provide a unified self- 
consistent description of the sequential and cotunncling 
transport regimes in the case of a molecular switch in the 
"classical" regime V > ojo and eV > Tiluq. This regime 
allows for a systematic non-perturbative treatment for an 
arbitrary electron-ion coupling strength.— We determine 
the dynamics of the vibrational degree of freedom, the 
average current and current noise through the device, and 
the switching times as functions of transport and gate 
voltages. We also analyze the role of extrinsic dissipation. 



II. MODEL 

We consider the model for a molecular switch pro- 
posed in Ref. Hllll5 . The molecule is modeled as a single 
electronic level d strongly interacting with a vibrational 
mode, x. It is located between two leads, from which 
electrons can tunnel into the electronic level. The inter- 
action is provided by the force A (typically of electrostatic 
origin) acting on the molecule. The system is described 
by the Hamiltonian 



H = 



(e + Xx)d^d+^- 
2m 



E 

k . at 



(i) 



taicU+Scka), (2) 



where a is the lead index (L or R) and c and d are the 
electron annihilation operators for the leads and local 
orbital, respectively. We consider the model for spinless 
electrons for simplicity. (Inclusion of spin along with on- 
site Coulomb blockade should lead to qualitatively sim- 
ilar results.) The vibrational mode is characterized by 
the "bare" frequency ujq and the effective mass m. The 
displacement and coordinates are described by the canon- 
ically conjugate operators x and p. The coupling between 
the electronic level and the mode is characterized by the 
"polaron" energy E p = A 2 /(2mw 2 ) and the coupling to 
the leads by tunnel rate r Q = irv a t 2 r /h, where v a is the 
density of states in lead a. In Refs. [lllfl2l it has been 
shown that for strong enough coupling, E p /K ^> Fl+Th, 
the system can exhibit bi-stability, with one state cor- 
responding to empty resonant level and non-displaced 
mode x, and the other to occupied level and the mode dis- 
placed by the amount ~ A/(tocl> 2 ). In the previous work, 
Ref. IT2I current and current noise were determined in 
the regime of small transport voltage, \eV\ <C E p (where 
eV = /J,l — P-r) in the approximately "symmetric" situ- 
ation, eo ~ E p . In the present work, we generalize the 
previous results for current and current noise as well as 
determine the behavior of the switching rates between 
the metastable states for arbitrary transport and gate 
voltages. 

When electrons are driven out of equilibrium by an ap- 
plied transport voltage, the dynamics of the vibrational 
mode becomes very simple, even for strong coupling be- 
tween the mode and electrons. That is because when the 



characteristic timescale for electronic subsystem becomes 
shorter than oscillator frequency c^o, electrons appear to 
the mode as a "high-temperature," albeit position de- 
pendent and strongly coupled bath. Physically, for any 
position x, the electronic bath adjusts (almost!) instanta- 
neously, in a manner analogous to how electrons adjust to 
the instantaneous positions of ions in isolated molecules, 
as described by the Born-Oppenheimer approximation. 
Indeed, as in the standard Born-Oppenheimer approxi- 
mation in equilibrium bulk solids, one effect of the non- 
equilibrium fast electronic environment is the modifica- 
tion of the effective potential experienced by the mode; 
however, what is more, the electronic subsystem, by 
virtue of being open, also provides force noise (fluctu- 
ations) and the dissipation to the mode. Since the force 
acting on the mechanical mode is simply —An, where 
n = d*d is the occupation of the electronic mode, in 
order to obtain the average force and its fluctuation it 
is enough to calculate the average of n and it's fluctua- 
tion (charge noise) for a given static position x. When a 
weak time dependence of x(t) is included one finds that 
a correction to the average of n appears that is linear 
in dx/dt. This last term corresponds to the dissipation 
induced by the retardation of the electronic degrees of 
freedom, that do not respond immediately to a change of 
x (first non-adiabatic correction)^ It can also be traced 
to the "quantum" nature of the charge noise, i.e. a slight 
asymmetry between the charge noise at positive and neg- 
ative frequencies i 30 ! 31 ' 32 As a result, the dynamics of the 
mode x becomes essentially classical, described by the 
Langevin equation, 12 



rax. + A(x)x + mux^x = F(x) + 



(3) 



where the position-dependent force F, damping A, and 
the intensity of the white noise D, — D(x)S(t— 

t') are related to the electronic Green functions on the 
Keldysh contour as 



F(x) 



x 2 h 



d(jjGf r (uJ, x), 



(4) 



M x ) = -ir~ I dujG fr {uj,x)d UJ G r f(uj,x), (5) 
Ztt J 

D{x) = — — / dwGf r (oj,x)G r f(uj } x) . (6) 
2n J 

The zero temperature Green functions (for the forward- 
reverse Keldysh time path) are 



G f r (w,x) 

G rf (uJ,x) 



hT L Q(^ L - Tiuj) + hT R <d([i R - hid) 
2i — : ,(7) 



-2i 



(hu - e - Ax) 2 + h 2 T 2 
hTL<d(huj — (Xl) + ?iTrQ(Tiuj — fin) 
(Hcu - e„ - Ax) 2 + H 2 T 2 



Here r = Tr, + Rr. These expressions are valid also at 
finite but low temperatures such that fc^T < hT. [At 
higher temperature the step functions 0(e) have to be 
replaced by Fermi functions np(— e/fcsT).] Therefore, at 
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low temperatures, we obtain, 



F(x) 



A(x) 



X 

+ T R ^fam 
X 2 Th 3 



e — Xx ir 
" + 2 



HT 

i P-R — e — Xx 7r 
KT + 2 



^ I [(Mi - eo - A.t) 2 + ^r 2 ]^ 

Ifl 

[(^ - eo - Xxf + h 2 T 2 ] 2 



D(x) 



x 2 r L r 



L± R 



Trr 3 



tan 1 z + 



Z 2 + 1 / VR-tQ-X, 



(9) 

(10) 
(11) 



Note that the expression for the force is just F — —Xn(x), 
where n(x) is the occupancy of the d level for a fixed 
displacement x. The expression for D is given for /i£ > 
fiR, otherwise, the p,L and [Ir have to be interchanged. 



III. CURRENT AND NOISE FROM THE 
FOKKER-PLANCK DESCRIPTION 

From the Langevin Eq. ^ one can derive a Fokker- 
Plank equation for the probability P(x,p,t) that at a 
given time t the displacement and the momentum of the 
vibrational are x and p = mx, 



dtV = --d x V - F(x)d p V ■ 
m 



A(x) 



d P ( P P) 



D(x) 



(12) 

This Fokker-Plank equation can be used to study both 
the stationary properties of the system, as well as the 
time evolution from a given initial condition. 



A. Current 

Given our assumption about the separation between 
the slow ionic - vibrational - and fast electronic - tunnel- 
ing - timescales, the problem of evaluating the stationary 
current reduced to the evaluation of the quasistationary 
current averaged over the fast electronic times for a fixed 
position x and momentum p of the mode, with the con- 
sequent averaging over the stationary probability distri- 
bution, P(x,p). In our case, the quasistationary current 
through the molecule depends then only on the position 
x (for k B T < KT), 



1 { x ) = TT \ duT(u,x) 



with 



T(uj,x) 



4r L r, 



(w - e - Xx) 2 + V 2 



(13) 



(14) 



The expectation value current is then simply 

I(t) = dxdpV(x,p)I(x) . (15) 



Solving the stationary Eq. (fT2|) one can thus obtain the 
current voltage characteristics for the device. 



B. Current noise 

We are also interested in the current noise: 

S(to) = f dfe*-* (l(t)J(O) + J(0)I(i)\ , (16) 



where 1 = 1— and / is the current (quantum) oper- 
ator. Again, since in our problem we have a clear time- 
scale separation between the vibrational and electronic 
degrees of freedom, we can distinguish two contributions 
to the current noise. The first is quasistationary (for 
a given position x) shot noise which arises due to the 
discrete nature of the electron charge. It has the usual 
form for a device with a single channel and transparency 
T(x,oj)B, 



S s hot(u = 0,x) 



2e^ 

h 



div 
2^ 



T(u, x)[l - T(u>, x)] . 



(17) 

The only change due to the presence of the oscillator is 
the fact that it must be averaged over the position, in the 
same way as we have done for the average current above. 

The second more interesting type of noise is caused 
by the fluctuations of the position x. It occurs on a 
long time scale, and thus, at low frequencies, it can be 
much more important than the standard electronic shot 
noise. 28 When the typical electronic and mechanical fluc- 
tuation times are of the same order of magnitude one 
has to take into account the correlation between the two 
sources of fluctuations,— However, for our system the 
separation of the time scales makes these two noises ad- 
ditive and allows for their separate evaluation without 
regard for one another. 

To obtain the low frequency "mechanical" contribution 
to the noise one needs to consider the autocorrelator of 
the quasistationary current (|15p at different times. This 
requires knowledge of the time-dependent solution of the 
Fokker-Plank equation (TT2"]) . The evolution of the prob- 
ability can be rewritten in a more compact form as 



d t V = CP 



(18) 



where C is the Fokker-Planck operator, in this notation 
V is a vector (Pi) and £ is a matrix (Ai)- The index 
i = (x,p) represents all the stochastic variables in dis- 
crete notations. For instance, the current operator X is 
diagonal in the i variables [cf. Eq. (fTS")) ] so that the av- 
erage current can be written simply as 



(J) = ^liVQi = (wo,Zvq) 



(19) 
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where v n i and w n i are the right- and left-eigenvectors of 
C with eigenvalue X n (Cv n — X n v n and w\C — A„wjj). If 
the eigenvalues are not degenerate then one can always 
choose the normalization so that (w n ,v m ) = 8 n>m . The 
conservation of the probability implies that Ao = 0, and 
by definition vo is the stationary solution and woi = 1. 
The fluctuation operator for the current is X = X — (X) 
in terms of which we can define the current fluctuations: 



S{t>0) = Y,^U l j{t)X J v Oj 



(20) 



Here Uij(t) is the conditional evolution probability that 
the system evolves from the state j at time to the state 
j at time t. It must satisfy the evolution equation (|18p 
with the boundary condition Uij(0) = 5ij. By Laplace 
transform [U(s) = f °° U(t)e~ st dt with Re s > and 



U(t)=r+;™(ds/2m)U(s)e 



a > 0] we obtain 



(s - C)U(s) = U(t = 0) = 1 



(21) 



We can then calculate the noise spectrum by using the 
symmetry S(t) = S(—t), 

S(lu) = S(s = -iu) + 0+) + S(s = iuj + 0+) (22) 

where S(s) is the Laplace transform of S(t) and has the 
form 



v 0j 



(23) 



We thus obtain 



S(u) = -2^i, 



c 



c 2 



X jVQj . (24) 



IV. RELEVANT PARAMETER RANGE 

We assumed from the beginning that wo<cr. This en- 
sures that the electronic dynamics of the device is faster 
than the vibrational one. The only remaining relevant en- 
ergy scale is E p , which we have to compare to the other 
two parameters Tiloq and %T. If L ^ E p /h the switching 
effects are difficult to observe since the boundaries of the 
Coulomb diamonds arc blurred on a scale HT much larger 
than the energy scale of the vibrational motion. We thus 
will not investigate this limit, but shall concentrate on 
the opposite one of T <C E p /Ti. 

It is convenient at this point to rewrite the Fokker- 
Planck equation in dimensionless form by introducing the 
variables y — kx/X, r = tujo, q = pk/\ujQm. Eq. (fTJ 
becomes 



-qd v V 



Td q V 



Ad q (qP) 



2 q 



(25) 



with 

T(y) = -y - 1/2 
+7H tan - 
where ji — F,;/F. 



1l tan" 



x (v g -v/2-y 



x (Vg+v/2-y 



(26) 



Ay) 



^r 2 



1L 



+ 



1R 

[K-1>/2-tf) 2 +f2]3 



, + u/2- H 



7T r V Z 2 + 1 J Vg-v/2-y 

We have also introduced the bias and gate voltages 



(27) 



(28) 



(Ml + Mfl)/2 - e = 2v g E p (29) 



and the dimensionless system parameters F = 2HT/E p 
and Cj — 2Thoq/E v . 

We can now discuss the limit of interest ujq <C L -c 
E p /h. The fluctuating and dissipative parts of the 
Fokker-Planck equation (coefficients A and T>) are much 
smaller than the force term (JF) since they are propor- 
tional to Cj <^ 1. For Cj —> the force term remains finite, 
while A and T> vanish. One therefore expects that the 
evolution of the system can be further coarse-grained in 
time. The system evolves under the influence of T most 
of the time and thus conserves its effective energy defined 
by E eff(y, q) = U ef f(y, q) + q 2 , with 



Uef f (y) = - / dy'Hy') 



(30) 



The effect of the small terms A and T> is to produce a slow 
drift among the nearby constant-energy orbits. The sta- 
tionary solution should then be a function of E e ff(y,q) 
alone and it is possible to reduce the Fokker-Planck equa- 
tion to an energy differential equation that in presence of 
a single minimum has the analytical stationary solution 



Q(E) = Me 



= Kf J E (a(E')/fj(E'))dE' 



/(3(E)- 



(31) 



Here j\f is a normalization factor and Q(E,t) 
j dydqS(E — E e ff(y, q))V(y, q, r). The coefficients a and 
P are obtained by averaging a combination of A and 
V on the trajectories of given constant effective energy 
E e ff(y,q) = E , as discussed in detail in Ref. [25|: 
a = (V(y)/2 - A{y)q 2 ) E and (3 = (p 2 V(y)/2) E . Note 
that in Eq. (f3"Tj) Cj cancels out in the exponential. Thus 
the limit Cj — > is well defined for the stationary distri- 
bution of probability. Obviously in this limit the time 
to reach the stationary state diverges since it is linearly 
proportional to Cj. 
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FIG. 1: Regions in the v — v g plane of existence of the minima 
of Ueff for r — > 0. The letters A, B, and C stand for the 
presence of a minimum at y — — 1, — 7l, and 0, respectively. 
The plane is separated into three dashed regions according to 
which of the three extrema is the absolute minimum. 




FIG. 2: Regions of stability of the sequential tunneling so- 
lution for jl ~ 1/2 (left panel) and 7l = 0.1 (right panel), 
and f = 0.02, 0.04, 0.08, 0.16, and 0.30, (from red to blue). 
The region of sequential tunneling evolves from the small tri- 
angular shape in the top of the plot (for T small) to a large 
trapezoidal shape (for large F) that touches the v = axis. 
The regions to the left and the right of the sequential tunnel- 
ing are "blocked" in the or 1 occupation state, respectively. 



When the potential can be approximated by a 
quadratic function around a local minimum and the y 
dependence of the coefficients A and V can be neglected, 
the expression for the probability becomes 

Q{E) = Ne- E ' T * , (32) 

where T* = 2T>{y m ) / A{y m ) and y m is the position of the 
local minimum. 

Even if in the general case the stationary distribution 
is not determined in such a simple way it is instructive 
to study the structure of U e ff(y). This is particularly 
simple for T <C 1 since in this limit the force becomes 

Hv) = -V - 1L 0(v g + v/2 -y) - lR 6{v g - v/2 - y) . (33) 

It is then possible to show that the effective potential 
landscape can show up to three minima at the positions 
y = for v g < -v/2, y = -~/ L for -v/2 - -y L < v g < 
v/2 — 7^, and y = for v g > v/2 — 1. (For simplic- 
ity we consider only the v > case.) The minimum at 
y = — Jl is due to the sequential tunneling for which the 
average occupation of the dot is < 7l < 1 (the energy 
level lies in the bias window). The other two minima 
correspond instead to classically-blocked transport (thus 
co-tunneling is the dominant current mechanism), either 
in the n = or n = 1 state. There are regions where two 
or three minima are present at the same time. One can 
show that for — v/2 — 7l/2 < v g < v/2 + 1 — 7^/2 and 
v > 1/2 the sequential tunneling minimum at x = —"/l 
is the absolute minimum. In the rest of the plane either 
the blocked state 0, or the blocked state 1 are true min- 
ima, the separation line between the two joins the point 
v g = — 1/2, v = to the apex of the conducting region 
v g = -3/4 + 7jj/2 and v = 1/2. (cf. Figs. Hand [21) 
For finite value of T the stability diagram changes, the 
main difference is the increase of the region of sequential 
tunneling that extends towards the axis v = 0, as shown 
in Fig. d 



At low voltage and small f one of the two blocked 
states has the minimum energy. For ji = = 1/2 
and v g = —1/2 the effective temperature of these states 
vanishes linearly with the bias voltage. Thus for v — > 
these are the "cold" states. The effective temperature 
at the sequential tunneling minimum (x = —1/2) is 
T* = 7ri> 4 /2 4 /r, thus for small T this state is always 
"hot." Around v = 1/2 the hot sequential tunneling state 
becomes the U e tf minimum, and the system starts to 
fluctuate between the hot and cold states. The dimen- 
sionless current / = I /Ye in the cold state is very small 
~ Tv while in the sequential tunneling regime it is of 
the order one. The fluctuations between these two states 
produces large telegraph current noise, as discussed for 
small v in Ref. H3 . 

The fact that the effective noise temperature varies as 
a function of the position can lead to dramatic conse- 
quences. In the conventional equilibrium statistical me- 
chanics, according to the Gibbs distribution, the low- 
est energy state is the most probable one. However, if 
the noise temperature varies as a function of position, it 
may happen that the lowest energy state, if it experiences 
higher temperature, may be less likely than a higher en- 
ergy state that experiences lower temperature. We illus- 
trate this point in Fig.[3l which compares the naive effec- 
tive potential profile U e f / with the actual self-consistent 
probability distribution. 

We need to stress here, however, that we assume that 
the only environment that is experienced by the mechani- 
cal mode so far is the non-equilibrium electronic bath due 
to the attached leads. If the dominant environment were 
extrinsic (non-electronic), with a fixed temperature and 
the coupling strength, then the effective potential would 
indeed uniquely determine the probabilities of particular 
states. We will come back to this point in Section [Villi 

In order to discuss the behavior of the device in the 
full range of parameters here we resort to a numerical 
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FIG. 3: Effective potential U e ff(y) (red dashed) compared 
to u ef f = -\nV(y) (blue) for f = 0.08, j L = 0.1, G) = 
10~ s ,v g — and different values of v as indicated in the panes. 
The quantity u e f / plays the role of an effective potential if T* 
was constant. Note in particular the case v — 1.2 for which 
the absolute minimum of U e ff is not the absolute minimum 
of u e ff due to the fact that T* is much lower in the other 
minimum. 



solution of the Fokker-Plank equation from which we can 
determine both the current and the current noise of the 
device. In the following section we discuss the numerical 
results. 




FIG. 4: Current for F = 0.02, 0.04, 0.08, 0.16, and 0.30, 
from the lowest to the highest curve at low bias. The other 
parameters are u) = 10 -3 , jl = 1/2, and v g — —0.5. 
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V. NUMERICAL RESULTS FOR THE 
CURRENT AND ZERO-FREQUENCY NOISE 

The expressions (|T5|) . (|17l) and ([24)) can be used to cal- 
culate the current and the noise of the device. In general 
the analytical evaluation of these expressions is not possi- 
ble. Numerically, the solutions can be obtained by rewrit- 
ing Eq. (fl"2"|) on a discrete lattice (x,p) and replacing the 
derivatives with their finite differences approximations. 
If the equation is solved in a sufficiently large (e.g. rect- 
angular) region in the x-p plane, one can use vanishing 
boundary conditions, since the probability vanishes far 
from the origin. The matrix corresponding to the dis- 
cretized Fokker-Plank operator C is very sparse and the 
numerical solution is relatively easy for matrices of di- 
mensions up to 10 5 . The discretization step sizes kAx/ X 
and Apk/Xujom must be smaller than HT/Ep in order 
to have a good convergence. This practically limits our 
numerical procedure to values of HT / E p > 0.01. 

We begin by considering the symmetric case, "/l = 1/2. 
The current as a function of the voltage bias for different 
values of T is shown in Fig. 0J One can see that for 
f — > the current is suppressed for v < 1/2 and rises very 
rapidly for transport voltages exceeding the threshold, 
as expected from the qualitative arguments given above. 
Numerically is difficult to reduce T further, but we expect 
that for r — > a discontinuity should appear as found in 



FIG^ 5: Fano factor of the current noise in logarithmic scale 
for f = 0.02, 0.04, 0.08, 0.16, and 0.30, from the lowest to 
the highest curve at large bias. Co = 10 -3 , 7l = 1/2 and 

Vg = -0.5. 



the case when cotunnelling is negligible. 25 

In Fig. we plot on a log scale the Fano factor 
[F = S(lu = 0)/2ei] of the mechanically generated cur- 
rent noise (the standard shot noise contribution is much 
smaller). One can see that F reaches huge values of the 
order of 10 3 , while it is typically 1 for the purely elec- 
tronic devices. The maximum of the Fano factor appears 
slightly below the value of the voltage where there is a 
crossover from the the cold to the hot minima; we will 
see later that this corresponds to the value for which 
the switching rates between the two minima are nearly 
the same. Since the blocked minimum is colder than 
the sequential tunneling minimum, this crossover hap- 
pens before the hot minimum becomes a true minimum. 
Enhancement of noise in this device should serve as a 
strong indication of the presence of mechanical oscilla- 
tions. 

In Fig. [51 [3 [5] and [5] we show the behavior of the Fano 
factor in the plane v g — v for T — 0.08 and 7^ = 0.5 
or 0.1. Note that in the asymmetric case, Fig. [5] and [51 
there is a very sharp peak in the Fano factor if we increase 
the bias voltage at fixed gate voltage greater than zero. 
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FIG. 6: Symmetric case. Fano factor for the induced cur- 
rent noise as a function of v g and v for T = 0.08, = 0.5 
and u> = 10 -3 . 
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FIG. 9: Asymmetric case. Current and Fano factor for the 
mechanically induced current noise as a function of v g and v 
for f = 0.08, 7l = 0.1 and Co = 10" 3 . 

This structure appears at the threshold of the sequential 
tunnelling conducting region. 

VI. SWITCHING RATE 




0.5- 




FIG. 7: Symmetric case. Current and Fano factor for the 
mechanically induced current noise as a function of v g and v 
for f = 0.08, 7l = 0.5 and Q = 10" 3 . 
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FIG. 8: Asymmetric case. Fano factor for the mechanically 

; a; 

7l = 0.1 and Co = 10" 3 . 



In the previous sections we have studied the current 
and the current noise. These quantities are the most 
readily accessible in transport measurements; however, 
it is interesting also to investigate what is the typical 
switching time t s between the two minima. This quantity 
can give an indication if the telegraph noise could be 
detected directly as a slow switching between discrete 
values of the average current. For this to happen the 
switching time must be very long-at least comparable to 
the average current measurement time (typically, in the 
experiment > 1 fis) . 

To find a reliable estimate of r s we need to know the 
typical time necessary for the system to jump from a 
local minimum of the effective potential Eq. ([50]) to a 
neighboring one. This concept is well defined since the 
diffusion and damping term of the Fokker-Plack equation 
are very small and the time evolution of the system on 
a short time scale is controlled by the drift term. Let 
us denote the value of the effective potential at the local 
maximum separating the two minima of interest as E max . 
The region f2 on the y-q plane around the minimum de- 
fined by E mm < E e ff(y,q) < E max can be considered 
as the trapping region. If the system is at time at the 
position (y, q) inside we can estimate the average time 
to reach the boundary of Q (dQ,) by solving the equation: 



C T 



-1 



(34) 



with (absorbing) vanishing boundary conditions on dfl. 35 
Here r stands for the function r(y, q). Since we are inter- 
ested on the average time to leave the region we average 
the escape time with the quasi-stationary distribution 
function. The vanishing boundary conditions introduce 
a sink thus there is no zero eigenvalue for the C operator 
with vanishing boundary conditions on dft. We can nev- 
ertheless always identify the eigenvalue with the smallest 
real part and call it Aq: Cvq — XqVq. We thus obtain 



(r, vp) 

(Mo) 



l 



(35) 
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FIG. 10: Symmetric case. j L = 1/2, f = 0.08, v g = -1/2 
and cli = IO -3 . Switching time between the two minima: red 
full line for the blocked transport minima (y = or y = — 1), 
and blue dashed line for the sequential tunneling minimum. 
In the inset: the current in each minimum (same notation of 
main plot) the average current (black full line) and the current 
noise (magenta dashed line). 
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FIG. 11: Asymmetric case. -y L = 0.1, f = 0.08, v g = -1/2, 
and ui — 10 -3 . Same notations as as in Fig. 1101 



The inverse of the lowest eigenvalue gives the average 
switching time; this is not surprising since the time evo- 
lution of the eigenstate vo is e~ tx °. It decays exponen- 
tially on a time scale — 1/Ao due to the escape at the 
boundaries of the region f2. 

We implemented numerically the calculation by solv- 
ing the Fokker-Planck equation in the energy-angle coor- 
dinates. If (y , 0) is a minimum of the effective potential 
with energy E m i n , we rewrite the Fokker-Planck equa- 
tion in terms of the variables E(x,q) — q 2 /2 + U e ff(x) 
and 9(x, q) = arctan((j/ (x — x )). In this way the bound- 
ary conditions read V(E = E max ,9) — for all values 
of 9. The results are shown in Figs. [TU] and QT] for the 
symmetric and asymmetric case, respectively. 

Let us begin by discussing the symmetric case of Fig. 
1101 For small bias voltage only two minima are present, 
they are perfectly symmetric and they correspond to two 
"blocked" (classically-forbidden) current state with n = 



or 1. The switching time is very long, and the system 
switches between two blocked states, each with very small 
cotunncling currents. Since the cotunneling currents for 
both minima in the symmetric state are the same, there 
is no telegraph noise for small v. As it can be seen from 
the value of result for the noise, the current fluctuations 
are nevertheless high, and the reason is that to jump 
from one minimum to the other the system has to pass 
through a series of states for which current flows through 
the device is significant. Moreover the slow fluctuations 
of the distribution function inside each minimum are im- 
portant for the noise as discussed in the following Sec- 
tion I VIII The fact that the jumping times are so long 
may actually hinder the observation of the jumps in a 
real experiment with finite measurement time. In a real 
device then the noise could be smaller in that case. In- 
creasing the voltage to v « 0.28, the sequential tunneling 
minimum at x = — .5 appears and a true telegraph noise 
start to be present. We see very clearly this in the escape 
times, which are no longer symmetric (we plot the y = — 1 
and y = —1/2 minima escape times, the y — minimum 
has the same behavior of the y = — 1 minimum), and 
the average current at the minima also changes abruptly. 
Even if the noise has a strong maximum near v — 0.28 
there is not a dramatic increase at the appearance of 
the minimum. The presence of the cotunneling smoothes 
the transition also for the noise that has its maximum 
before the sequential minimum appears. The switching 
time changes by 6 orders of magnitude in a very small 
range of bias voltage. Above v w 0.53 only the sequential 
tunneling minimum survives. 

We consider now the asymmetric case of Fig. QTJ It 
is clear that the evolution of the escape times is very 
different from the symmetric case. In particular we con- 
sider the strongly asymmetric case of 7l = 0.1. In this 
case the sequential tunneling minimum merges with the 
blocked n = minimum, leading to a two minima land- 
scape of the potential. The consequence is that there is 
no abrupt appearance of a new minimum for some value 
of the bias voltage; rather, the two minima are always 
present at the same time till v « 0.4. At low voltage the 
potential landscape is nearly symmetrical, both minima 
are cold, but for the sequential tunneling one is charac- 
terized by a slightly higher T* and thus its escape time is 
shorter (dashed line in Fig. ITT)) . Increasing the voltage, 
the height of the potential barrier for the blocked state 
reduces, thus reducing the escape time. At some point 
(in the case of Fig. QT]for v ~ 0.18) the escape time from 
the cold state become shorter than the escape time of the 
hot one, since the the temperature has to be compared 
with the barrier, and at this point the barrier height is 
smaller in the cold state. Near the crossing region the 
noise shows a maximum, due to the fact that the system 
spends nearly half of his time in each of the two minima, 
with different average current. Tuning v one can thus 
cross from a region where the system is trapped in one 
of the two minima, to a region where it jumps on a rela- 
tively long time scale from one minimum to the other. If 
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FIG. 12: Average Current (black continuous line), Fano fac- 
tor (magenta dot-dashed line), current in the two minima 
(same notation as Fig. [10} for f = 0.08, j L = 0.1, w = 10" 3 
at v g — 0. In the inset the escape times from the two minima 
as a function of the bias voltage v. 



the switching time scale becomes of the order of the re- 
sponse time of the measuring apparatus it is in principle 
possible to observe directly the fluctuation between the 
two values of the current. 

This is even more pronounced if we follow the evolution 
of the current at v g = 0. As can be seen in the contour 
plot of the Fano factors (cfr. Fig. [9J, in this way we will 
cross a very sharp peak of the Fano factor. The results 
are shown in Fig. [T2J At low voltage only a single nearly 
blocked state is present (x = — 1 and n = 1). For v w 0.8 
a new minimum appears at x « ~1l — —0.1 that is 
for the moment at higher energy and with a very small 
barrier. The current associated to this minimum is much 
higher than the other, and the system starts to switch 
between the two states. The switching is very slow thus 
the noise is high. Very rapidly as a function of v the new 
local minimum becomes the true minimum, and then the 
other minimum disappears. 



VII. FREQUENCY-DEPENDENCE OF THE 
CURRENT NOISE 

Eq. derived above can be applied to study not 

only the zero frequency noise, S(uj = 0), and the Fano 
factor, as we did in Section [Vl but also the current noise 
at an arbitrary frequency. In this section we numerically 
evaluate S(uj) and provide a qualitative explanation for 
the observed trends. As we mentioned before, the shot 
noise contribution to the noise can be neglected as far as 
the frequency considered is much smaller than T. From 
the numerical calculations we find that the frequency de- 
pendence is characterized by a single frequency scale, and 
approximately is Lorentzian peaked at u) = 0. This can 
be seen in the inset of Fig. [T3l where we show S(u) as a 
function of lo /luq on a logarithmic scale for several values 
of the bias voltage v. One can parameterize each curve 
by a single number, that we choose as the frequency lo c 



at which S(lo c ) — S(0)/2. It is instructive to compare 
the time scale 1/lo c with the energy dissipation and the 
switching timescales in various regimes. 



At low voltages, since switching between the 
metastable minima is exponentially slow, we anticipate 
that the low frequency (lo < loq) current fluctuations will 
be determined by the energy fluctuations within the sin- 
gle well in which the molecule spends most of its time. 
For a simple harmonic oscillator, the corresponding time 
scale is given by the inverse damping coefficient. For 
small energy fluctuations, the current changes with en- 
ergy linearly. Thus, current fluctuations will track the en- 
ergy fluctuations, i.e. will be Lorentzian with the width 
given by A/m. To check this we plotted in Fig. [T3l the 
value of rmoo / A(x) evaluated at the minimum of the po- 
tential (dotted line) . There is a reasonable agreement for 
low voltage but, as expected, not for large voltages. The 
reason is that at large v the system becomes hot, and the 
energy dependence of A cannot be neglected. To address 
this issue, we calculated the average of A(x) with the dis- 
tribution function V(x) obtained by solving numerically 
the stationary problem. The result using thus obtained 
A is shown as dashed line on the figure. We find that it 
agrees very well with the lo c extracted from the numer- 
ical calculation of S(ui), both at high and low voltages. 
Note that at high voltage the energy dependence of A is 
crucial to understand the frequency response of the noise. 
The effective temperature changes the average of A, and 
hence lu c by nearly three orders of magnitude. 



In the intermediate transport voltage regime, 1 < v < 
1.3, the system switches between the two wells frequently. 
Therefore, we naturally expect that the timescale for the 
current noise should depend on the switching rate be- 
tween the wells. If each of the wells would corresponds 
to a fixed value of current the resulting noise would be a 
telegraph, with the Lorentzian lineshape and width given 
by the sum of the switching rates. However, in each well 
as a function of energy current is not fixed. In fact, the 
current increases gradually in the "blocked" well as the 
energy the approaches the top of the barrier reaching the 
value / ~ r near the top of the barrier. On the other 
hand, in the well where transport is sequential, current 
remains approximately / ~ T for any energy. There- 
fore, one can naturally expect deviations from the simple 
telegraph behavior. Indeed, we find that the timescale 
1 /lo c tracks the escape time from the "blocked" well (blue 
dot-dashed line in Fig. [T3")) . which is the longer escape 
rate, and the fast escape from the "hot" sequential well 
does not matter. We therefore conclude that the noise is 
governed by the energy (and thus current) fluctuations 
within the cold (more probable) well, which also occur 
on the timescale comparable to the escape rate from it. 
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FIG. 13: Inset: Frequency dependence of the current noise 
for several values of the bias voltage. From this data we es- 
tracted u) c as the frequency at which S{u) c ) = S(0)/2. Main 
plot: comparison of loo/uj c red full line, the escape time u>qt 
blue dot-dashed line, the friction coefficient uJom/A at the 
minimum dot line, and averaged light dashed line. The pa- 
rameters are the same as Fig. 1121 
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FIG. 14: Voltage dependence of the current and Fano factor 
(inset) for different values of the extrinsic dissipation: fj = 0, 
10 -3 and 0.1. The temperature of the external bath is 0.01 in 
our dimensional units, the other parameters are the same as 
Fig. 1121 One can see that the current jump becomes sharper 
for stronger coupling to the environment. At the same time 
the Fano factor becomes sharper, thus a strong noise region 
survives, but becomes very narrow when the external bath 
dominates. 



VIII. ROLE OF EXTRINSIC 
ENVIRONMENTAL DISSIPATION 

As we discussed above most of the effects we found are 
due to the non-equilibrium dynamics of the oscillator. In 
order to improve our understanding of this fact, and to 
probe robustness of the results to external perturbations 
we consider the influence of extrinsic dissipation on the 
system. This can be easily included in the model since 
the coupling to an external bath implies only additional 
dissipation and fluctuation on top of the intrinsic ones. 
We assume that the system is damped due to the cou- 
pling to an external bath at equilibrium at the temper- 
ature Tfc. The fluctuation and dissipation coming from 
this coupling satisfy the fluctuation dissipation theorem. 
Thus the presence of the extrinsic damping induces the 
following change in the variables A and D defined in Eqs. 
([TP]) and CP: A -> A + rj and D -> D + k B T b ri/2. We 
present the numerical results for the dimcnsionlcss pa- 
rameters fj = rj/muiQ and Xb = 2ksT/Ep. The numerical 
procedure remains unchanged. 

We show in Fig. Q3]the behavior of the current and the 
noise for the same parameters of Fig. [P^lbut at X& = 0.01 
and for different values of the external dissipation. The 
main feature that can be clearly seen is the sharpening of 
the step for the current. The external damping reduces 
the position fluctuations of the oscillator thus reducing 
its ability to escape from the blocked regions of the pa- 
rameters' space. On the other side if the oscillator is in a 
conducting region the probability that it can fluctuate to 
regions of blocked transport is smaller, thus the current 
is increased in the conducting regions and reduced in the 
blocked regions, increasing the steepness of the step. For 
the same reason the region of large noise is reduced. We 
find that the value of the Fano factor remains actually 
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FIG. 15: Effective potential U e ff (red dashed) compared 
to — ln(V(x)) for different values of the extrinsic dissipation: 
fj = 10" 5 , 10~ 3 ' 5 , 10" 2 , and 10" ' 5 , from the lowest to the 
highest curve. The temperature of the external bath is 0.01 
in our dimensional units, the other parameters are the same 
as Fig. [12] and the curves are shifted and multiplied by a 
constant factor for clarity. 



very large, but only in a very narrow range of bias volt- 
ages. Increasing the coupling to the external bath reduces 
this windows and thus finally rule out the possibility of 
observe it at all. 

A second interesting quantity to study is the distribu- 
tion function V(x). If the coupling to the environment 

dominates we expect that V(x) = const e~ UB ff( x '< Tb to 
verify this fact we compare U{x) = — \nV{x) and U e ff(x) 
in Fig. [TS] We find that for small coupling first U(x) de- 
viates even more from the form of U e ff. the minimum 
in the cold regions deeps (left minimum in the figure). 
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The reason is that the increase of the damping is more 
effective in the cold region where both of damping and 
fluctuation are small. In the hot region (right minimum 
in the figure) the intrinsic fluctuation and dissipation is 
very large and for small external damping there is no 
noticeable effect. Increasing the coupling to the environ- 
ment also the hot minimum is cooled and the shape of 
U becomes similar to that of U e ff shown dashed in the 
plot. This shows how relevant the non-equilibrium dis- 
tribution of the position is for the determination of the 
transport properties of the device. 

IX. CONCLUSIONS 

In this work our goal was to provide a unified descrip- 
tion of the transport properties of the strongly coupled 
non-equilibrium electron-ion system mimicking a molec- 
ular device, in a broad range of parameters. Our results 
are based on a controlled theoretical approach, which 
only assumes that the vibrational frequency is the low- 
est energy scale in the problem. In this regime, the vi- 
brational mode experiences the effect of the electronic 
environment as a non-linear bath that has three inter- 
related manifestations: (1) Modification of the effective 
potential, including formation of up to two additional 
minima, (2) position-dependent force noise that drives 
the vibrational mode, and finally, (3) position-dependent 
dissipation. We have self-consistently included the effect 
of tunneling electrons on the dynamics of the vibrational 
mode, and the inverse effect of the vibrational mode on 
the electron transport. This enabled us to obtain the 



average transport characteristic of the "device," i.e. the 
dependence of the current on the transport and gate volt- 
ages, as well as address the problem of current noise and 
mechanical switching between the metastable states. The 
agreement between the switching dynamics and the fre- 
quency dependence of the current noise determined inde- 
pendently, enabled us to construct a comprehensive but 
simple understanding of the combined electron-ion dy- 
namics in different transport regimes. In particular, the 
enhancement of current noise may serve as an indicator of 
generation of mechanical motion, and its magnitude and 
frequency dependence provide information on the regime 
the molecular switching device is in and the values of 
relevant parameters. 
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